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Abstract 

We introduce a modification of the Fast Marching Algorithm, which solves the anisotropic 
eikonal equation associated to an arbitrary continuous Riemannian metric on a two or 
three dimensional box domain. The algorithm has a logarithmic complexity in the maxi- 
mum anisotropy ratio k(M) of the Riemannian metric M., which allows to handle extreme 
anisotropies for a reduced numerical cost. We establish that the output of the algorithm 
converges towards the viscosity solution of continuous problem, as the discretization step 
tends to zero. The algorithm is based on the computation at each grid point z of a reduced 
basis of the lattice 7L d , with respect to the symmetric positive definite matrix A4(z) encoding 
the desired anisotropy at this point. 

Introduction 

The eikonal equation, and its generalization the Hamilton-Jacobi equation, is a Partial Differen- 
tial Equation (PDE) which describes an elementary front propagation model: the speed of the 
front depends only on the front position and orientation. This PDE is encountered in numerous 
applications, such as motion planning control problems |19j . modeling of bio- medical phenom- 
ena [IT], and image analysis [15] . It was also recently used in the context of medical image 
analysis [1] for extracting vessels in two dimensional projections or three dimensional scans of 
the human body, and for performing virtual endoscopies. This application requires to solve a 
highly anisotropic generalized eikonal equation with a high resolution on a cartesian grid, at a 
computational cost compatible with user interaction. It is one of the key motivations of this 
paper. 

This paper is devoted to the construction and the study of a new algorithm, Fast Marching 
using Lattice Basis Reduction (FM-LBR), designed to solve the anisotropic eikonal equation 
associated to a given Riemannian metric M, and which can handle large or even extreme 
anisotropies. The domain must be of dimension two or three, unless the metric has a special 
structure, see Point hi at the end of this introduction, and discretized on a cartesian grid. 
FM-LBR, as its name indicates, is a variant of the classical Fast Marching algorithm \19[ 122] . 
an efficient method for solving the eikonal equation when the metric is isotropic (proportional 
at each point to the identity matrix). Lattice Basis Reduction [T3] is a concept from discrete 
mathematics, discussed in the first section of this paper, and used in the FM-LBR to produce 
local stencils for the discretization of the eikonal equation. Lattice Basis Reduction is involved 
here for the first time in the numerical analysis of a PDE to the knowledge of the author. We 
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describe the algorithm itself in the second section, where we also analyze its computational cost, 
and establish its consistency: the discrete approximations produced by this algorithm converge 
towards the viscosity solution of the continuous eikonal equation as the discretization step tends 
to zero. We present some numerical experiments in the third section, which confirm the small 
numerical cost of the algorithm and show an accuracy competitive in general, an remarkable in 
test cases related to the envisioned medical application. 

We consider a positive integer d, the dimension, which is fixed throughout this paper. Let us 
mention that our algorithm only applies in dimension d G {2, 3}, unless the problem of interest 



has a special structure, see Point [m] at the end of this introduction. We denote by VL the periodic 
unit box 

Q := (lR/7Z) d . (1) 

The domain Q, coincides with the unit box [—1/2, l/2] d equipped with periodic boundary con- 
ditions. We denote by z + v G fi, the offset of a point z G Vt by a vector v G H d . The periodicity 
assumption is not essential, as shown in the numerical experiments and discussed in Point [n] at 
the end of this introduction, but it simplifies the description and the proof of convergence of our 
algorithm. We consider a fixed Riemannian metric 

M G C°(Q,S+) 

where <Sj~ denotes the set of d x d symmetric positive definite matrices. For each u G H d and 
each M e Sj we denote \\u\\m '■= v 7 u T Mu. Our objective is to compute (an approximation 



of) the viscosity solution D : O, — > H + , see [IT] and Definition 2.9, of the anisotropic eikonal 
equation 

||VD(.z) = 1 f° r almost every z G ft \ {0}, . . 

D(0) = 0. ' [} 

The only element of general theory used in this paper is the uniqueness of the viscosity solution 
to the above problem. Another point of view on this problem is provided by the following 
characterization for each z G f2, the quantity D(z) is the length of the shortest path 7 

joining z to the origin. More precisely the solution of the eikonal equation is 

D(z) = D(z,0), (3) 

where D(-, •) denotes the Riemannian distance, defined for all x, y G f2 by 

D(z,y) := inf{length( 7 ); 7 GC7 1 ([0,l],O),7(0)=x, 7 (l) = y}, (4) 

length(7) := / h' (t)\\ M{l{t)) dt. (5) 
J 

We define the anisotropy ratio k(M) of a matrix M G S~t, and the maximum anisotropy 
k(M) of the Riemannian metric A4, as follows (denoting by || • || the standard euclidean norm) 

k(M):= max J^, k(M) := maxK(7W(z)). (6) 
[|u[|=|[u|[=l ||i>||m z ^ 

Note that k(M) = y/WMWWM' 1 ]]. The test cases presented in ^3] cover both moderate anisotropy 
k(M) ~ 5, and "extreme" anisotropy k(M) ~ 100, which is relevant in applications to image 
segmentation and structure extraction [3]. 

The periodic box Q is discretized on a cartesian grid Q n of step size 1/n 

n n ■- {0, 1/n, 2/n, • • • , (n - l)/n} d C O. 
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Figure 1: Some classical stencils used in the discretization of two dimensional (left) or three 
dimensional (right) eikonal equations, associated to a Riemannian metric A4. A Dijkstra inspired 
method, using one of these stencils at a point z G Q n , will be consistent if either Ai{z) is diag onal, 
or k(M{z)) < kq where Kq = 1, 1 + v2, 1, (V3+ l)/2 (from left to right). See Proposition 1.2 



Denoting by N := n d = #(Q n ) the number of grid points, our algorithm has complexity 

0(N]nN + NIxlk(M)), (7) 

where the constant behind the O(-) notation only depends on the dimension d. In practical 
applications one has N > 10 000 and k(A4) < 100, hence the first term dominates in (|7]), and 
the complexity of our algorithm is only 0{N In N). 

Let z G O and let 1/ C O be a compact neighborhood of z, not containing the origin in its 
interior. It follows from the interpretation of D as a path-length distance ([3j that 



D(z) = min B(z,z')+B(z'). 

z'aav 



(8) 



Our algorithm, as well as most alternative solvers of the eikonal equation, discretizes this PDE 
using an approximation of the right hand side of Q: the Hopf-Lax update operator [191 E]) 
earliest references in [51 E]. Equation Q is then reinterpreted as a fixed point problem. More 
precisely, consider a fixed n > 1 and assume that a small simply connected neighborhood V n (z) 
of each z £ £l n has been constructed under the form of a simplicial mesh (a triangulation if 
d = 2, a tetrahedral mesh if d = 3, see Figure [I]), with all its vertices in H n . 

We denote R + := R + U {+oo}, equipped with the topology of a compact segment, and we 
adopt the convention x oo = 0. For any discrete map d : fl n — > R + , the Hopf-Lax update 
operator is defined by 



An(d,z) 



mm 

z'edv 



z ~ z \\M(z) 



+ l V d(z') 



(9) 



where z E £l n is arbitrary, has neighborhood V := V n (z), and Iy denotes the piecewise linear 
interpolation operator on this simplicial mesh. The difference z' — z is well defined in ^ because 
z and z' both belong to the small and simply connected set V C Q, the stencil of z, which can 
be identified to a subset of TR d (this convention will not be recalled systematically in the rest of 
the paper). The Hopf-Lax update operator A n approximates the distance D(z, z') in Q by the 
local norm \\z' — ^||x( 2 ), for z' G dV, and the value D(z') by linear interpolation of d on the 
mesh V. Note that A n (d, z) depends on d(z') for all vertices z' £ (!„ fl dV, on the intersection 
of the discretization grid with the stencil boundary. 

The discrete approximation d ra : £l n — > H+, of the continuous solution D of the eikonal 
equation, is obtained as the solution of the iV-dimensional fixed point problem 



A n (d n ,z) = 
d„(0) = 0. 



d n (z) for all z G O n \ {0}, 



(10) 
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Numerical solvers of the eikonal equation differ by (i) the construction of the stencils V n (z), 
z G Q, n , and (ii) the approach used to solve the system ( |10[ ) which is inspired by the algorithms 
of Bellmann-Ford or of Dijkstra used in graph theory. The algorithm presented in this paper, 
FM-LBR, belongs to the category of Dijkstra inspired algorithms with static stencils, and among 
these is the first one to guarantee a uniform upper bound on the stencil cardinality. 



Bellman- Ford inspired algorithms. The discrete fixed point problem (10) is solved 
via Gauss-Seidel iteration: the replacement rule d(zk) <— A n (d, Zk) is applied for k = 
0,1,2,... to a mutable map 5 : Q n — > Hi, until a convergence criterion is met. In the 
fast sweeping methods, see [21] and references therein, the sequence of points (zk)k>o 
enumerates repeatedly the lines and the columns of Q n . Alternatively this sequence is 
obtained via a priority queue in the Adaptive Gauss-Seidel Iteration (AGSI) of Bornemann 
and Rasch [§J. The stencil V(z) of a point z £ Q n is usually the offset by z of a fixed 
stencil V given at the origin, such as those illustrated on Figure [TJ 

Fast sweeping methods have 0(\(A4)N) complexity when the metric Ai is isotropic (pro- 
portional to the identity at each point), but this result does not extend to anisotropic 
Riemannian metrics, see [26] for the proof and the expression of A(A4). The AGSI has 
complexity 0(n(A4)N 1+ d^ for arbitrary anisotropic Riemannian metrics, where is 
a non explicit constant which depends on global geometrical features of the metric j5]. The 
AGSI is a popular, simple, and quite efficient method, which is included for comparison in 
our numerical tests. 



Dijkstra inspired algorithms. The system ( 10 ) is solved in a single pass, non-iteratively, 



using an ordering of £l n determined at run-time. This is possible provided the Hopf-Lax 



update operator satisfies the so-called "causality property", see Lemma 2.3, which can be 
ensured if the stencil V(z) of each z £ Q n satisfies some geometrical properties depending 
on Ai(z), see Definition |1.1| The different Dijkstra inspired methods are characterized by 
the construction of the stencils V(z), which depend on Ai(z), in contrast with Bellman- 
Ford inspired methods which are characterized by the choice of the sequence (zk)k>o- Solv- 



ing the system (10) with a Dijkstra inspired algorithm has complexity C(/x(A / ()A r In N), 
where n(M) is an upper bound for the cardinality of the stencils (the number of simplices 
they are built of). 

In the Ordered Upwind Method (OUM) of Sethian and Vladimirsky |19l [2~5*] . the sten- 
cils are constructed at run-time; their cardinality is bounded by (D(K(A4) d ) and drops to 
0(n{M.) d ~ 1 ) as N — > oo. In contrast, the stencils are constructed during a preprocessing 
step and then static in the Monotone Acceptance Ordered Upwind Method (MAOUM) 
of Alton and Mitchell [2]; their cardinality is bounded by 0(n(M) d ). The FM-LBR in- 
troduced in the present work uses an approach similar to the MAOUM, except that the 
cardinality of the stencils is 0(1), fully independent of the Riemannian metric Ai. The 
complexity estimates are thus 0(n(M) d N\n N) for the OUM and the MAOUM (asymp- 
totically O^MY^NhiN) for the OUM), and 0(N\nN + Nhm(M)) for our approach, 
FM-LBR, where the second term in the complexity accounts for the stencil construction. 

The above mentioned algorithms are consistent for the anisotropic eikonal equation associ- 
ated to an arbitrary continuous Riemannian metrics Ai : Q — > in the sense that the discrete 
output d n of the algorithm converges to the viscosity solution D of the continuous problem as 
n — > oo. Some more specialized variants of the fast marching algorithm are only consistent 
for a restricted set of metrics, but can be executed nonetheless with an arbitrary anisotropic 
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metric Ai (In that case the discrete system (10) may not be solved, and the numerical results 



are variable, see £j3j). For instance the original fast marching algorithm [22\ is consistent if A4(z) 
is proportional to the identity matrix for each z £ £1, and more generally if Ai(z) is a diagonal 
matrix. In addition to these cases of isotropy and axis-aligned anisotropy, some variants [7] are 
also consistent if k(M) < Ko, where hq is a given bound, see Figure [l]for some examples in two 
and three dimensions. Our numerical experiments include for comparison one of these methods: 
Fast Marching using the 8 point stencil (FM-8, center left stencil on Figure [TJ, which is popular 
in applications [1] thanks to its short computation times and despite the lack of convergence 
guarantee for arbitrary metrics. Depending on the implementation I19 | [22], involving either 
a sorted list or a bucket sort, these methods have complexity O(NlnN) or 0(T(Ai)N), where 



?(M) := /max||X(z)||max||^(z') 



In the applications for which our method is intended, one typically has ln(iV) < k(A4) < 
T(Ai) <C N, in such way that the complexity ([T]) of the proposed method is comparable to 
O(NlnN) and smaller than 0(T(M)N). In summary the algorithm proposed in this paper has 
a complexity O (N \n N + N In k(M)) which is significantly less than general methods such as the 
AGSI, the OUM and the MAOUM, which are consistent for arbitrary Riemannian metrics, and 
no more than the complexity of more specialized methods, such as the original fast marching 
algorithm. 

Our complexity analysis guarantees short and predictable run-times, which is desirable in 
applications involving user interaction, e.g. image processing. It does not answer however the 
question of accuracy. 

Remark 1 (Accuracy and stencil construction). The above compared solvers of the eikonal 
equation, the AGSI, the OUM, the MAOUM and our method the FM-LBR, involve different 
stencils and thus a different Hopf-Lax update operator A n . The discrete fixed point problem 
therefore depends on the method, and so does its solution d n , which impacts the accuracy of 
the scheme. In the special case of a constant Riemannian metric, the theoretical error analysis 
presented in [12], shows that the anisotropy ratio k{M) does not impact the accuracy of FM- 
LBR, in an average sense over grid orientations. In contrast one expects the numerical error of 
alternative methods to grow proportionally to k(M), as discussed in 119^ and suggested in 
by numerical experiments for the AGSI. This error analysis does not extend to general Rieman- 
nian metrics, and we do not claim that the proposed algorithm outperforms its alternatives in 
terms of accuracy in all applications. The numerical experiments presented in section |3] show 
nevertheless that FM-LBR is competitive in this regard. 

Let us emphasize that the efficiency of the proposed algorithm, the FM-LBR, comes at the 
price of its specialization. We review below its limitations: 

i (Finsler metrics) FM-LBR only applies to the anisotropic eikonal equation associated to 
a Riemannian metric. The underlying Riemannian structure plays an important role in 
our approach, and our algorithm therefore cannot handle more general Hamilton-Jacobi 
equations, contrary to the AGSI [5], the OUM [19] and the MAOUM [2] mentioned above. 
These approaches allow Finsler metrics in addition to Riemannian metrics, in other words the 
local euclidean norm || • ||x( 2 ) m ay be replaced for each z £ £1 with an arbitrary asymmetric 
norm | • \ z 6 C°(lR d ,]R+). 

Constructing static stencils suitable for Dijkstra inspired algorithms, in the sense that they 
guarantee the causality property as in Lemma |2.3| (the stencil is said to be causal), is an 
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active subject of research in the case of Finsler metrics. A characterization was obtained 
in [23] , and a construction proposed in [21 [13] . It was also observed in [TBI COS E] that the 
canonical stencil (Figure [TJ left) is causal for a large variety of axis-aligned Finsler metrics, 
in such way that the original fast marching algorithm can be applied. 

ii (Domain discretization) FM-LBR requires a domain discretized on a cartesian grid. Our 
algorithm therefore does not apply to domains provided under the form of general un- 
structured meshes, a more difficult setting which has attracted an important research effort 
E [HJ, [TU1 [19] . At the price of a higher technicality of the proof, the specific periodic domain 
(llj) could however be replaced with a more general smooth bounded domain C TR d , still 
discretized on a grid. In that case the PDE (12]) is replaced with 



where the boundary data / needs to satisfy [TT] the Lipschitz regularity condition \f(z) — 
f(z')\ < length(7) for any path 7 G C 1 ([0,l],O) joining two points z,z' G dQ. Note that 
using the FM-LBR in this context will require the extension of the boundary data to a ghost 
layer covering dtl, since this algorithm uses large stencils, of euclidean diameter 0(n(M)/n) 
instead of e.g. C(l/n) for the AGSI. 

iii (Dimension) FM-LBR applies to arbitrary continuous Riemannian metrics on a two or three 
dimensional domain (an extension to four dimensional domains is presented in |12j). The 
algorithm can be extended to higher dimension if the Riemannian metric Ai has a specific 
diagonal block structure, with blocks of size 1, 2 or 3, see §2. Such block diagonal structures 
are not uncommon in the context of medical imaging, see [3J. They are inherited from the 
cartesian product structure of the fast marching domain: f2 = £Iq x tlx, where £lo is a physical 
domain of dimension < 3, and Q\ is an abstract parameter domain of dimension < 2. 

We introduce and study in ^ljthe notion of M- reduced mesh, where M G Sj[ is a symmetric 
positive definite matrix. The construction of M-reduced meshes of bounded cardinality, using 
Lattice Basis Reduction, is the main originality of this paper and the key of the small complexity 
of our algorithm the FM-LBR, which uses them as stencils. Following a more classical approach, 
we describe in ^2] the FM-LBR as a variant of the Fast Marching algorithm using this specific 
stencil construction, and we establish the related convergence result. We finally we present some 
numerical experiments in §[3] 

Remark 2 (Computing distances on a surface). Consider a smooth surface S C IR? (or more 
generally a smooth embedded manifold S C IR d ), equipped with Riemannian metric induced by 
the euclidean metric on IR? . In order to compute distances on S, a first possibility is to solve 
an isotropic eikonal equation on S, by applying the OUM on a triangulated mesh of S fWlj , or 
the original fast marching algorithm if this triangulation only contains acute triangles. A second 
possibility is to use one or several local charts ip : O — > S and to solve an anisotropic eikonal 
equation on Q; grid discretisations are allowed since n C M 2 , hence the FM-LBR can be applied, 
as well as numerous alternatives. In practice the data usually dictates the choice between these 
approaches, since changing from one representation to the other (triangulation or charts) is a 
non-trivial and computationally intensive procedure. 




for almost every z G 17 
for all z G dVl. 



(11) 
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1 Reduced meshes and bases 



We introduce in this section the notion of M-reduced meshes, where M G ST is a given symmet- 
ric positive definite matrix. Anticipating on §[2j we motivate their study by their use as static 
stencils for Dijkstra inspired eikonal solvers, such as the FM-LBR. Consider an Riemannian 
metric M. G C°(Q, Sj~), and assume that for each z in the grid Q n a A4(z) -reduced mesh T(z) 
has been constructed. Denote by V n (z) the stencil obtained by rescaling and offsetting the mesh 
T(z): with obvious notations 

V n (z) = z + -T(z). (12) 
n 

We show in ^2] that the discrete fixed point system (10), obtained with these stencils, can 
be solved in a single pass, and we establish a convergence result towards the solution of the 



continuous eikonal equation, in Theorem 2.2 



A simplex T C H a is the convex hull of d+ 1 points vq, ■ ■ ■ , Vd G not lying on a common 
hyperplane, which are called the vertices of T. A mesh is a finite collection T of simplices which 
satisfy the following conformity condition: for all S,T G T the intersection S n T is the convex 
hull of the common vertices of S and T. 

Definition 1.1. A M-reduced mesh is a mesh T which satisfies the following properties. 

(a) The union of the elements of T is a neighborhood of the origin 

(b) The vertices of each simplex T G T lie on the lattice 7L d , and T has volume l/dl. 

(c) For each T G T , one of the vertices ofT is the origin 0, and the others denoted by vi, ■ ■ ■ ,Vd 
satisfy for all 1 < i < j < d 

vjMvj > 0. (13) 
Heuristically Point (a) of Definition 1 1 . 1 1 ensures that the numerical information is propagated 



in all directions in the discrete fixed point system (10). Point ( b) i mplies that the non-zero 
vertices of any simplex T G T form a basis of 7Z d , see Definition 1.4 a property later used to 
show that the information does not "fly over" a subset of the grid f2 n . Point (c) is an acuteness 
condition related to the celebrated Causality property, see [19], which is a prerequisite for the 
Fast Marching algorithm and other one-pass solvers of eikonal or Hamilton- Jacobi equations. 

The next proposition gives a simple criterion to show that a given mesh T is M-reduced for 
all M G Sj" of sufficiently small anisotropy k(M). 

Proposition 1.2. Let T be a d-dimensional mesh which satisfies the requirements (a) and (b) 
of Definition Let 



K 01 '■= \ ] +7 !3-l ' where y(T) := min - ^ , (14) 
y 1 - j(T) T,(u,v) \\u\\\\v\\ 

and where the minimum in j(7~) is taken a non-zero vertices u,v of a common simplex T G T ■ 
The mesh T is M-reduced for any M G S"j~ such that k(M) < k(T). 

Proof. Let u,v be two non-zero vertices of a common simplex T G 7~, and let M G ST. Let 
u' := u/\\u\\ and let v' := v/\\v\\. By construction we have 

\\u' + v'f = 2(1 + u' T v') > 2(1 + 7 (70), IK - v'f = 2(1 - u'V) < 2(1 - 7 (T)). 
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Let us assume for contradiction that u T Mv < 0, which implies that \\vf + v'\\m < \\ u ' — v'\\m- 
Observing that 



' „,'I|2 



K ( M) 2 = HMiiiiM- 1 !! > rr^ ! U ', +V X > ^t3> 

v ' - \\u'-v'\\ 2 \\u' + v'\\ 2 M 1-7(T) 

we obtain that k(M) > k(T), which concludes the proof of this proposition. □ 

The meshes illustrated on Figure [I] are easily checked to satisfy requirements (a) and (b) of 
Definition Requirement (c) clearly holds for diagonal matrices. Evaluating (14), we obtain 
k(T) = 1,1 + a/2, 1) (V3 + l)/2 (from left to right), hence (c) also holds for matrices M G S j~ 
of anisotropy k(M) within this bound. These meshes are therefore .M(z)-reduced under the 
conditions described in the caption of Figure [TJ hence the resulting system of equations (10) can 
be solved in one pass by a Dijkstra inspired method. 

Remark 1.3 (Block diagonal matrices). Let di,d,2 be positive integers, let Mi G S^, and let 
M2 G S^ 2 . Let T\ be a M\-reduced mesh, and let 75 0, Mi-reduced mesh. Let d := d\ + di and 
let T be the d-dimensional mesh defined as follows: for any T\ G 71 of vertices = Uq, • • • u^, 
and any T2 G 75 of vertices = Uq, ■ ■ ■ , u^ 2 , the d-dimensional simplex T of vertices 

(0,0),(u}, ()),••• ,«,0),(0, «?),••• ,(0,<) 

belongs to T '. Then one easily checks that T is a M -reduced mesh, where M G denotes the 
matrix of diagonal blocks Mi and Mi- 



We introduce in subsection £1.1 the algebraic notion of M- reduced basis of 7L d , where M G 
and 1 < d < 4. We show that the collection of vertices of a M-reduced mesh contains a 



M-reduced basis, a property later used in the analysis of our algorithm. We next give in §1.2| an 
explicit construction of a M-reduced mesh of bounded cardinality, for each M G St, d G {2, 3}. 

1.1 Bases of the lattice TL d 

The results of this section fall in the framework of low-dimensional lattice basis reduction. We 
refer to |14j and references therein for an introduction to this rich theory, from which we use 



only one result: Theorem 1.5 stated below. 

We denote by u{7L + • • • + u^TL the sub-lattice of TL d generated by ui, ■ ■ ■ , u& G 7L d : 

u\7L + • • • + U]~7L := {u\Z\ + • • • + u^z^; z\, ■ ■ ■ , z^ G 7L}. 

If k = then the above sum equals {0} by convention. 

Definition 1.4. A basis of 7L d is a d-plet (ui, • • • , Ud) of elements of 7L d such that 

|det(m,--- ,u d )\ = 1. (15) 

Assume that 1 < d < 4. A M-reduced basis of7L d , where M G , is a basis (u\, ■ ■ ■ ,Ud) of7L d 
which satisfies for all 1 < k < d 

u k G argmin{||z||M; z G 7L d \ { Ul 7L + ■■■ + u k -t7L)}. (16) 
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Figure 2: The unit sphere {u; \\u\\m 
M-reduced mesh constructed as in Proposition 
ranging from 1 to 15, and eigenvector (cos(37r/8), sin(37r/8)) associated to the small eigenvalue. 



1}, a M-r educed basis (u,v), and the boundary of a 
for some M E of anisotropy ratio k(M) 



1.9 



For instance the canonical basis (ei, • ■ ■ , e^) of H d , is also a basis of 7L d . If M E St is a 
diagonal matrix of coefficients (Ai, • • • , Ad), and if < < • • • < ^a(d) f° r some permutation 
a, then (e CT (i), • • • e a (d)) is a M-reduced basis of 7L d . See Figures [2J |3| [4] and [5] for some examples 
of M-reduced bases associated to non-diagonal matrices M. 

A M-reduced basis of 7L d contains by definition small vectors with respect to the norm || • ||m : 
the smallest linearly independent ones with integer coordinates. As a result, the construction 
of M-reduced meshes based on M-reduced bases, presented in the next subsection, will produce 
small meshes, again with respect to the norm || • ||j^f. As a result, the discretization of the 



FM-LBR uses small stencils V n (z) (12) at grid points z E f2 n , in the sense of the local norm 
II • llxl 2 )- Small stencils often the key of good accuracy in the discretization of PDEs. This 
observation is the starting point of |12j . an error analysis of the FM-LBR, in the case of constant 
Riemannian metrics. 

The above definition of a M-reduced basis of 7L d is adequate only up to dimension 4, since 
in dimension d > 5 there exists matrices M E St such that no basis of 7L d satisfies the relations 



(16), see [H] (these relations state that ||wi||m equals the i-th Minkowski's minimum Aj(M)). 
The proper generalization of Definition 11.41 in dimension d > 5, is Minkowski's reduction 



Theorem 1.5 (Nguyen, Stelhe, 2009). There exists an algorithm which, given a matrix M E 
as input, 1 < d < 4, produces a M-reduced basis of7L d and has the numerical cost O (l+ln k(M)) . 



Proof. The proof is contained in |14j , and we only point out here the precise reference within the 
paper and the slight differences in notations. The algorithm described in [13] takes as input a 
basis (61, • • ■ , bd) (here: the canonical basis of TR d ) of a lattice L (here: 7L d ), and its Gram matrix 
with respect to some scalar product (here: the Gram matrix is M). The algorithm outputs a 
greedy reduced basis of the lattice L, a notion which coincides with Minkowski's reduction if 
d < 4 (Lemma 4.3.2 in [H]), which itself coincides with Definition 1.4 if d < 4. 



The main loop of the iterative algorithm is executed at most the following number of times 
(Theorem 6.0.5 in [Tl|): 



O ( 1 + In max 

Ki<d 



%\\M 



In min IliilU/ 
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hence 0(1 + In ||M|| 2 - In ||M _1 1|" 2 ) = 0(1 + 1iik(M)) times in our setting. The complexity of 
each of these iterations is dominated by a closest vector search, described in Theorem 5.0.4 in 
|14| . which consists of the inversion of a k x k Gram matrix, where 1 < k < d— 1, and a 0(1) 
exhaustive search. In terms of elementary operations (+, — , x, /) among reals, each iteration of 
this algorithm thus has cost 0(1), and the overall cost is the number of iterations 0(l+ln k{M)). 

Note that an important part of the discussion in [14J is devoted to the special case where 
the vectors (61, - • • ,bd) have large integer coefficients, the Gram matrix is computed with re- 
spect to the standard euclidean scalar product, and the complexity of an elementary operation 
(+, — , x,/) among integers is not 0(1) but depends on the size of these integers. This more 
subtle notion of complexity, named bit complexity, is not relevant in our setting. □ 



In dimension d = 2, the algorithm mentioned un Theorem 1.5 mimicks the search for the 
greatest common divisor of two integers, and is often referred to as Gauss's algorithm [T3] . This 
algorithm uses only a pair (u, v) of (mutable) variables in 2Z 2 , initialized as the canonical basis 
of R 2 . The pair (it, i>) becomes a M-reduced basis at the end of the following loop, which takes 
at most 0(hiK(M)) iterations. Round denotes rounding to a closest integer. 

Do (it, v) ::= (i>, it — Round(« T Mv/||u||j^) v), while \\u\\m > IMIm- 

Proposition 1.6. Assume that 1 < d < 4. Let M £ S"f and let (ui, ■ ■ ■ ,Ud) be a M-reduced 
basis of 7L d . Then for all 1 < i < d 



\\m\\ < k(M), 
\ui\\m < k(M)\\ui\\ M - 



(17) 
(18) 



For any integer combination z of the elements of the basis distinct from m, in other words 
z = a\U\ + • • • aj_iUi_i + ai+iUj+i + ctdUd, where a±, • • • , oti—x, a «+i> ' • ■ ,<J(i£S, one has 



2\ujMz\ < 



\?\\ 2 
\ z \\m- 



(19) 



Proof. We consider a fixed 1 < i < d, and we claim that 



II M _1 1 



1 

~ 2 < \\ui\\m < 



u. 



\\m < \\M\ 



(20) 



The left inequality follows from 1 < ||ui|| < \\M 2 1| ||iti||Af- The central inequality follows f rom 
the fact that u\ minimizes the norm || • \\m among all elements of 7L d \ {0}, see Definition 1.4l 
Denoting by (e\, ■ ■ ■ , e<f) the canonical basis of TR d , we observe comparing dimensions that there 
exists 1 < j < d such that e,- ^ ui7Z + ••• + ••• Ui-\7L. Therefore ||itj||M < ll e jlUf < ll-^ll 5 using 



Definition 1.4 which establishes (20), hence also (18). We obtain (17) combining (20) with the 

w w w 



observation 



< \\Ui \\m\\M~ 



We next turn to the proof of ( 19 ), and for that purpose we remark that Ui + z £ u\7L + • • • + 
Ui-\7L. Indeed otherwise itj would be a linear combination of u\, ■ • ■ , itj+i, • • • , Ud, which 



contradicts (15). Definition 1.4 thus implies that 



IKHm < IK + z\\\i = IKIlIf + ZujMz + ||z||!f 



which implies that —2ujMz < \\z\ 



M- 



We obtain likewise 2ujMz < 



proof of (19). 



which concludes the 

□ 
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We say that point z is a vertex of a mesh T, if it is a vertex of one of the simplices T G 7~. 
We introduce a distance d x on the collection Sj[ of symmetric positive definite matrices, which 
is defined as follows: for all M, N G 

d x (M, N) := sup | In \\u\\m — In IMIjvI • 

The distance d x allows to compare the norms of vectors multiplicatively, justifying the x sub- 
script, in contrast with the classical operator norm which is tailored for additive comparisons. 
Indeed denoting a := d x (M,N) and j3 := \\Mz —N^\\, one has for all u G M 2 such that [|u|| = 1 

e~ a < \\u\\m /\\u\\n < e a , and — f3 < \\u\\m — \W\\n < ft- 

The next lemma establishes a lower bound on the \\-\\m norm of points outside of a N- reduced 
mesh, when the matrices M, N G Sf are close enough. 

Lemma 1.7. Assume that 1 < d < 4. Let M, iV G Si". Let iti, • • ■ , 6e an arbitrary M -reduced 
frasis o/ TL , and let T be a N -reduced mesh. Consider a point z G TL which is not a vertex of 
T '. Then there exists 1 < I < d such that 

z£ Ul 7L + --- + ui7L and \\z\\ 2 M e Mx(M ' N) > \\ Ul \\ 2 M + \\ Ul \\ 2 M . 

Proof. Since the union of the elements of T is a neighborhood of the origin, there exists a simplex 
T G T and a real A > such that Xz G T. Denoting by v%, ■ ■ ■ ,Vd the non-zero vertices of T, 
there exists therefore non-negative reals «i, • • • , G IR + such that z = ai«i + • • • + a^Vd- Since 
| det(i>i, • • • , Vd)\ = d\\T\ = 1, the coefficient aj is an integer for all 1 < i < d. 

Up to reordering the vertices vi, ■ • • , u^, we may assume that fi, • • • , ffc are positive, and 
Vk+i,"- ,Vd are zero, for some 1 < k < d. We denote by I the smallest integer such that 
Vi, ■ ■ ■ , Vk G u\ZL + • • • + u{ZL. By additivity, z G U\2L + • • • + u{ZL. By assumption there exists 
1 < i < k such that Vi u\7L + • • • + ui-\7L, and thus ||i>i||m > ||^||m; other vertices satisfy 
H^illAf > 11^1 ||m; 1 < i < k, since they are non-zero and have integer coordinates. We thus 
obtain 

p 4d x (M,N) || _ II 2 > 2d x (M,iV) || ~||2 
e ll^llAf — M ^ M jv 

a 2 ||fi|| 2 v + 2 otiajvjNvj 

\<i<k l<i<j<k 



> e^xWTV) a *\\ Vl f N 



Ki<k 



> 



Y ^inim 



Ki<fc 



. 2<i<fc 



We have a\ + • • • + a\ > 2 since z is not a vertex of T, which concludes the proof. □ 

The following corollary shows that the vertices of a iV-reduced mesh contain a M-reduced 
basis, whenever M and N are close enough. This property will be used in the convergence 



analysis of the FM-LBR, Lemmas 2.6 and 2.7, to show the connectivity of the graph defined by 
the local stencils in the discrete domain. It is is illustrated on Figures [2] and [3| the basis vectors 
u, v are always among the vertices of the consecutive mesh, and conversely. 
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Figure 3: The unit sphere {u; \\u\\m 
a M-reduced mesh constructed as in Proposition 



1}, a M-reduced basis (u, v), and the boundary of 
for some M G 5^" of anisotropy ratio 



1.9 



k(M) = 6, and eigenvector (cos(0), sin(0)), 9 G [vr/4, 7r/2], associated to the small eigenvalue. 



Corollary 1.8. Assume that 1 < d < 4. Let M, N G Sj~ be such that 

d x (M, N) < ln(l + k(M)- 2 )/4. 



(21) 



Let (tii, • • • , Ud) be a M-reduced basis of 7L d , and let T be a N -reduced mesh. Then m, ■ ■ ■ ,Ud 
and —Ui, ■ ■ ■ , —Ud are vertices ofT. 

Proof. We consider 1 < I < d, and we assume for contradiction that ui (or —ui) is not a vertex 



of T. It follows from (21) and (18) that 



| N |2 fe 4d x (Af,A0 < ii^u^ + K ( M )- 2 \\ Ul \\ 2 M < \\ui\\ 2 M + Umlll,. 

The previous lemma thus implies ui £ u\7L + • • • + u^ZL for some 1 < k < I. This contradicts 
Definition 1.4, and concludes the proof. □ 



1.2 Explicit construction of M-reduced meshes 

This subsection is devoted to the explicit construction of a M-reduced mesh of bounded cardi- 

This 



1.10 



nality for any M G S^, where d = 2 in Proposition 1.9, and d = 3 in Proposition 
construction uses as a starting point a M-reduced basis of 7L d . 

Let A be a d x d invertible matrix. For any simplex T we denote A(T) := {Az; z G T}, and 
for any mesh T we denote A(J~) := {A(T); T G T}. 

Proposition 1.9. Let M G and let (u\,U2) be a M-reduced basis of 7L 2 . Let u := u\ and 

let v := £U2, where e G { — 1, 1} is chosen in such way that u T Mv > 0. Consider the collection 
T of triangles of vertices and one of the following pairs 

{u, v}, {v, v — u}, {v — u, —u} 

or their opposites (the opposite of the pair {a, b} is {—a, —b} ). Then T is a M-reduced mesh. 

Proof. We denote by (a),(b),(c) the corresponding points of the Definition |1.1| of M-reduced 
meshes. 

Consider a matrix A such that Au = (1, 0) and Av = (0, 1). Then the mesh A(T) does not 
depend on u and v, and the reunion of its elements is easily checked to be a neighborhood of 
the origin. This immediately implies (a). 
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Figure 4: A M-reduced mesh for particular a matrix M £ S£ , such that k(M) = 10 (left). 



Illustration for the proof of Proposition 1.9 (center left). Connectivity defined by ( |23| ) (center 
right) and by (p5|) (right), in Proposition 
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The equality of the determinants 

det(u, v) = det(v, v — u) = det(v — u, —u) 

shows that all the elements of T have volume 1/2. Since their vertices clearly lie on this 
establishes (b). 



We have by construction u^Mv > 0. Furthermore, applying (19) to the vectors u\ and 
z = U2 (resp. U2 and z = u\) we obtain 2\u T Mv\ < \\v\\ 2 M and 2\u T Mv\ < [|u[|^f, hence 

v T M(v -u) = \\vf M - u T Mv > 0, (v - u) T M(-u) = \\uf M - u T Mv > 0, 

which establishes point (c) and concludes the proof. 

For a visual illustration of point (c) we suggest the reader to look at Figure [4] (center left), 
which displays the image of the mesh T by a linear transformation P such that Pu = (1, 0) and 
P T P = M. For any two vectors a, b G H 2 , we have a? Mb > if and only if the images of a 
and b by P form an acute angle. The blue region corresponds to all possible values of v which 
satisfy the constraints 2\v T Mu\ < \\u\\ 2 M < \\v\\ 2 M . □ 

Proposition 1.10. Let M S and let [u\,U2,u^) be a M-reduced basis ofH? . We distinguish 
two cases depending on the parity of the number of non-negative scalar products among ujMu2, 
u^Mu^, ujMu\. 

• Odd parity. We denote u := e±ui, v := £2U2, w := £3^3, where £i,£2j £ 3 £ { — 1 , 1} are 
chosen in such way that 

u T Mv > 0, u T Mw > 0, v T Mw < 0. (22) 

Consider the collection T of tetrahedra of vertices and one of the following triplets 

{w,u,w + v} {w,w + v,w + v — u} {w,w + v — u,w — u} 

{w,w-u,-v} {w,-v,u-v} {w,u-v,u} , 23 - 

{v,w + v,u} {v, w + v — u, w + v} {v , v — u, w + v — u} 

{—u,—v,w — u} {—u,w — u,w + v — u} {—u,w + v — u,v — u} 

or their opposites (the opposite of the triplet {a,b,c} is {—a,—b,—c}). Then T is a M- 
reduced mesh. 
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Even parity. We denote u := E\u ax , u := E2U a2 , w := £3^0-3, where £1,62, £3 G {—1, 1} and 
the permutation a of {1,2,3} are chosen in such way that 

v T Mw > u T Mv > u T Mw > 0. (24) 

Consider the collection T of tetrahedra of vertices and one of the following triplets 



{u, v,w} {w,v,w — u} {w — u,v,v — u} 

{w — u,v — u,—u} {w — u, —u,w — v} {w — v,—u,—v} 

{w — v, —v, u — v} {w — v,u — v,w — v + u} {w — V + u, u — v , u} 

{u,w,w — v + u} {w,w — v,w — v + u} {w,w — u,w — v} 



(25) 



or their opposites (the opposite of the triplet {a,b,c} is {— a, — b, — c}. Then T is a Ad- 
reduced mesh. 

Proof. We denote by (a),(b),(c) the corresponding points of the Definition |1.1| of M-reduced 
meshes. 

In each case of the parity, if A is a matrix such that Au = (1,0,0), Av = (0,1,0) and 
Aw = (0, 0, 1), then the mesh A(T) is independent of u, v, w. The union of the elements of this 
constant mesh A(T) is easily seen to be a neighborhood of the origin, see Figure[4j which implies 
(a). 



The determinant of each triplet of vectors appearing in (23) and (25), is equal to det(u, v, w). 



Since | det(u,v,w)\ = 1, see (15), the volume of all the constructed simplices is 1/6. Furthermore 



their vertices obviously lie on 2Z 3 , which establishes (b). 

We next need to establish (c): the non-negativeness of the scalar products between two 
common vertices of any tetrahedron T £ T ■ To avoid notational clutter we denote in this proof 

(u, v) := u^Mv. 

We first remark that the positivity of the scalar products (u, v) , (v,v — u),(v — u, —u) (associated 
to edges which lie in the plane generated by (u,v)) follows as in Proposition |1.9| from the 
inequalities (u,v) > 0, 2\(u,v)\ < \\u\\\j and 2\(u,v)\ < \\v\\j^. For the other scalar products we 
need to distinguish between the two cases of the parity. 



Odd parity. It follows from (22) that 



(u,v) > 0, (u,w) > 0, (-v,w) > 0, 
(u, w + v) > 0, (w, u — v) > 0, (w — u, —v) > 0. 



Applying (19) to z = u we obtain 2\(u,w)\ < hence (— u, 



w — u) 
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\ 2 M ~ (u,w) > 0. 



Likewise 

{v, w + v) > 0, (w, w + v) > 0, (w, w — u) > 
Applying (19) to z = u we obtain 2\(u,v)\ < and 2\(u,w)\ < \\u\\ 2 M , which implies 



(— u, w + v — u) 



\ M — (u,v) — (u,w) > 0. Likewise 

(v, w + v — u) > 0, (w, w + v — u) > 0. 



Applying ( 19 ) to z = v — u we obtain 2\(w,v — u)\ < \\v — u||jL, hence (v — u, w + v — u 



\v — u\ 



M 



+ (w, v — u) > 0. Likewise 



(w — u,w + v — u) > 0, (w + v , w + v — u) > 0. 
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w=(2,l,l) 





Figure 5: The unit sphere {u; \\u\\m = 1}> a M -reduced basis (u, v, w), and a M-reduced mesh 
constructed as in Proposition 1.10, for the symmetric matrix M G Sg of eigenvalues 5 2 ,5 2 ,1 



(anisotropy ratio k(M) = 5) and eigenvector (4, 2, 3) associated to the small eigenvalue. 



Even parity. From the inequality (24) we obtain 



(v, w) > (u, v) > (u, w) > 0, (v, w — u) > 0, (—u, w — v) > 0. 

Applying (19) to z = u we obtain 2\(u, v)\ < \\u\l\j, and therefore (—u,w — u) > 0. Likewise 

(— v, to — v) > 0, (u;, to — u) > 0, (to, w — v) > 0. 
(u,w — v + u) > 0, (w, to — i> + it) > 0. 



Applying (19) to z = u we obtain 2|(u, u)| < [|«[||f and 2\(u,w)\ < ||«[||f, which implies 



(v — u,w — u) = (v, w) + (||mm|| 2 — (u, v) — (u, w)) > 0. Likewise 

(w — u, w — v) > 0, (w — v , u — v) > 0. 



Applying (19) to z 



v we obtain | (w, u — v)\ < \\u 



Likewise (w — v , w — v + u) > 0. 



\ M , hence (u — v, w — v + u) > 0. 

□ 



2 The algorithm 

We present in this section our modified fast marching algorithm, Fast Marching using Lattice 
Basis Reduction (FM-LBR). We estimate its complexity, and we establish its consistency: the 
convergence of the discrete approximations towards the solution of the continuous problem. Our 
presentation of this section is fairly classical, and we invoke (variants of) arguments previously 
seen in the literature on numerical schemes for the eikonal equation |22} [T9l [5]. A detailed 
description and proof is nevertheless necessary, since the original constructions of §1 do not 
exactly fit in earlier framework. The FM-LBR is split in two steps: Preprocessing (stencil 



construction), followed by Execution (Dijkstra-like resolution of (10)). The integer n > 2 is 
fixed in this subsection. 

After the Preprocessing of the FM-LBR comes its proper Execution. This part of the algo- 
rithm involves a boolean table b : £l n — > {trial , accepted} . Given such a table, we say that a 
grid point z E fi is accepted if b(z) = accepted, and that z is a trial point otherwise. 

Definition 2.1. Consider a map d : Vt n — > and a grid point x G We introduce a 

variant of the Hopf-Lax update A n (d,x) which involves two additional variables: a boolean table 
b : O n —> {trial, accepted} , and another grid point y G O n . 

A„(d, x; b, y) := min \\x' - x\\ M i x) + Iy d(x). (26) 
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FM-LBR: Preprocessing. 

Input. An integer n > 2. The values Ai(z), z G Q n , of a Riemannian metric Ai G C°(£l, S^). 



1. Produce a M (x)-reduced mesh 7~(x) for each x G fi n , using Proposition 1.9 in dimension 



2, Proposition 1.10| in dimension 3, or by combining these constructions with Remark 1.3 
in higher dimension if the metric has a block diagonal structure. 

2. Construct the stencils V n (x) = x + ^T(x), x G Sl„,, by offsetting and rescaling T(x). 

3. Compute the reversed stencils, defined by V n (y) := {x G Q n \{y}; y is a vertex of V^a;)}. 



VFe denoted by T the collection of vertices, edges and facets of the boundary of the stencil V := 
V n (x), which (i) contain only accepted grid points, and (ii) contain the point y. As before Iy 
denotes piecewise linear interpolation on the mesh V . The difference x' — x makes sense because 
x and x' belong to the small and simply connected stencil V of x. 
Equivalently, one has 

A n (d, x; b,y) := min < 

k,(ai),(vi) 

where the minimum is taken among all 1 < k < d, all a\ 1 ■ ■ ■ , at G M and all v\, - ■ ■ ,Vk G 7L d , 
such that 

• «j + • • • + afc = 1, and en > for all 1 < i < k. 

• Vx, ■ ■ ■ ,Vk are non-zero vertices of a common simplex T G T(x). 

• b(x + v\/n) = ■ ■ ■ = b(x + Vk/n) = accepted. 

• x + vi/n = y. 

The main originality of the FM-LBR lies in the stencil construction, in the Preprocessing 
step. The second part of this algorithm, Execution, is similar to other Dijkstra inspired solvers 
of the eikonal equation with static stencils, such as [21 122j . 

FM-LBR: Execution. 

4. Initialize to trial a (mutable) boolean table b : O n — > {trial , accepted} . 

5. Initialize to +oo another (mutable) table d : Q n — > H+, except for d(0) = 0. 

6. While there remains a trial point in f2 n 

(a) Denote by y the trial point which minimizes d, and set b(y) ::= accepted. 

(b) For all x G V n (y), set d(x) ::= min{d(x), A n (d, x; b,y)}. 

Output. The map d : Cl n — > IR + . 

The rest this section is devoted to the proof of the following theorem. Points (i) and (iii) 
also hold for any alternative construction of M. (z)-reduced meshes T(z), in the Preprocessing 
step, in dimension d < 4, which diameters are uniformly bounded by a constant C = C(k{M)). 



E 

Ki<k 



Vi 



n 



M(-A 



E 

Ki<fc 



a,- d 



(-3 



(27) 
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Theorem 2.2. (i) The map d n — d obtained at the end of the FM-LBR satisfies the discrete 



fixed point system (10) 



(ii) The cost of the FM-LBR is 0(N In N+N In k(M)), where N := n d = if elementary 

operations among reals (+,—,x,/ and \f-) have unit complexity. 



(Hi) Denoting by D : 17 — > 1R + the solution of the anisotropic eikonal equation 



one has 



lim ( max | D(z) — d n (z)| I = 0. 

n->oo \zefln ) 

The next lemma constitutes the causality property, which is the key to the one pass resolution 



of the system (10) in Fast Marching algorithms, see Proposition 2.4 



Lemma 2.3 (Sethian and Vladimirsky 2000, [19]). Let M G £j~ and let vi 

linearly independent vectors such that vj Mvj > for all 1 < i < j < k. Con. 
and convex set 



Choose A = (Si 



S := {a = (ax, ■ ■ ■ , ajfc) G R%] a\ H h a k = 1} 

, 5k) G M, denote 



5 := min 



OtiWi 

Ki<k 



M 



Ki<fe 



(28) 



and assume that the minimizer a G H o/ i/tis problem has positive coefficients. 

Then 5 > 5i for each 1 < i < k. Furthermore denoting by M the matrix of entries Mij 
vjMvj, and defining 1 := (!,••• ,1) G M k , one has the relations 



1 

Ma 



\S1-A\ 



M- 1 ' 



a 



M 



(51 - A). 



(29) 
(30) 



Proof. For completeness, and due to notation changes, we give the proof of this lemma, which 
follows the steps of Property A.l in [19J. 



The problem (28) is the optimization, on the compact an convex set H, of a strictly convex 



functional. Hence there exists a unique minimizer a G H. The gradients of the maps a i— > 



a 



M 



+ a A and a i— >■ a 1 , respectively the minimized function and the constraint in ( 28 ) , are 



proportional at the minimizer a according to Lagrange's Theorem. Hence there exists A £ R, 
the Lagrange multiplier, such that 



Ma 



+ A = Al. 



Therefore 



IA1-AI 



\a\ 



M 



Ma 



Using again (31) we obtain 
5 = Hal 



M 



+ a T A = a T 



a. 



Ma 



1. 



(31) 



(32) 



M" 1 



a 



+ A = a T (Al) = A. 



M 
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Combining the last equation with (32) (resp. plj)) we obtain ( p9[ ) (resp. (30)). 

Since M has non negative coefficients, positive diagonal coefficients, and since we have as- 
sumed that a has positive coefficients, we obtain that the product Ma has positive coefficients, 
hence also 51 — A. Therefore 5 > Si for all 1 < i < k in view of (30), which concludes the 
proof. □ 

The next proposition, adapted from [22], shows that the Execution of the FM-LBR solves 
the discrete fixed point system (10). 

Proposition 2.4 (Tsitsiklis, 1995). We denote by E the collection of maps d : Q n — ^ JR+ 
satisfying d(0) = 0. The set E is partially ordered as follows: d < d if d(z) < d (z) for all 
z G Q n . We denote by d° G E the map equal to +oo on f] n \ {0}. 

We denote by A n : E — > E the operator defined by A n (d)(z) := A n (d, z) if z G f2„ \ {0}. 

1. The operator A n is increasing on E: A n (d) < A n (d') for all d, d' G E such that d < d'. 

2. The sequence A^(d°) converges to a fixed point d G E of A n , as k — > oo. 

Let N := n d . We denote by (yi)o<i<N the grid points consecutively selected at step 6 (a) of the 
FM-LBR, and by (d*)i<j<Ar the maps obtained at the end of each iteration of the while loop. 

3. For all0<i<j<N one has d* > d j > d. 
4- For all < i < N one has d J (yi) = d(y,). 

5. For all < i < j < N one has d{yi) < d(yj). 

Combining Points 3 and 4 we obtain d N = d. Recalling Point 2 and observing that a fixed point 
of A n solves (10), we obtain Point (i) of Theorem 2.2. 

Proof. Point 1, the monotonicity of A n on E, follows from the expression Q of the pointwise 
entries A n (d, z), d £ E, z G il n , as a minimum of monotone maps (linear forms with positive 
coefficients). 

Point 2. Any map d G E satisfies d < d°, hence A(d°) < d°. Using Point 1 and an immediate 
induction argument we obtain that the sequence A fc (d°) is decreasing, and thus converging to a 
limit d G E. The identity A(d) = d follows from the continuity of A n , which is clear in view of 
its expression ^ as a minimum of linear functions indexed by a compact set. 

Point 3. Consider a fixed < i < N, and denote by b{ the value of the boolean map at 
the beginning of the i-th iteration (thus bi(yj) = trial if and only if i < j). The maps d l and 
d* +1 may only differ on the set V n (yi), see step 6 (b) of the FM-LBR. We assume that d l > d, 
an induction hypothesis that is satisfied for i = 0, and observe that d % > d i+1 by construction. 
Furthermore, for x G V n {yi) we have either d t+1 {x) = d L {x), hence d t+1 {x) > d(x) by induction, 
or 

d i+1 {x) = A n (d\x; b h yi) > A n (d,x; b hVi ) > A n (d,x) = d(x). 

Thus d l > d %+ > d, for all < i < N, which concludes the proof of Point 3. 

Following the steps of Lemma 3.2 in [22], we prove Points 4 and 5, simultaneously, by 
induction on i. Let y be the trial point (6j(y) = trial, equivalently y = yj for some j > i) which 
minimizes d. If y G dQ* = {0}, then y = 0, i = 0, d(0) = d°(0) = 0, and there is nothing to 
prove. Otherwise we have 

d(y) = A n (d,y) = \\y - ^ OLrVa{r)\\M{y) + ^2 a ^(y CT(r) ), (33) 

Kr<k Kr<k 
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for some 1 < k < d, positive barycentric coefficients (a r )^ =1 , and grid points (y a ( r ))r=i- These 
points satisfy d(y (7 ( r )) < d(y), due to the causality property, hence bi(y a ^) = accepted. Thus 
a(r) < i, for 1 < r < k, and therefore d* _1 (y .( r - ) ) = d(y c i r \) by induction hypothesis on Point 
4. Recalling the Hopf-Lax updates performed at step 6 (b) of the FM-LBR, in the previous 



iterations of the while loop, we find that the right hand side of (33) is an upper bound for d*(y), 
hence d l (y) < d(y). Since d* > d (Point 3), it follows that y minimizes d* among trial points, 
similarly to yi. Finally d(y) > d l (y) = d t (yi) > d(yi) > d(y), using Point 3 and the optimality 
of y for the last two inequalities. Thus d(yi) = d t (yi), and this is the minimum of d among the 
trial points {y,-; j >i}. This establishes Points 4 and 5 at rank i. □ 



We next turn to Point (ii) of Theorem 2.2 : the complexity estimate of the FM-LBR. Consider 



a point z G Q, n . Obtaining a A^(z)-reduced basis has cost 0(1 + In k(M. (z))), see Theorem 1.5 



and from this point constructing a Ai (z)-reduced mesh T(z) using Proposition 1.9 or 1.10 has 
cost 0(1). Step 1 of FM-LBR Preprocessing thus has total cost 0(N + N\uk(z)). Steps 2 and 
3 only cost 0{N). 



The FM-LBR Execution involves the evaluation of the Hopf-Lax update operator (27). Each 



mesh T(z), z G Q n , has 0(1) vertices by assumption, hence evaluating this operator amounts to 



solve 0(1) optimization problems of the form (28). Their solution is obtained as the root 5 of 



a univariate quadratic equation (29), hence the complexity is 0(1). Globally in the Execution, 



such evaluations happen at most the following number of times: 

Y J #(V n {z))= ^2n(z) = 0(N), 

where n(z) denotes the number of non-zero vertices of T{z). In our case n{z) = 0(1). In 
addition, the FM-LBR Execution requires to maintain a list of the elements of ft -1 (trial), 
sorted by increasing value of 5. Elementary insertions and deletions in this list have cost 0(ln N) 
(if this sorted list is implemented numerically using a binary heap), and occur each time the 
Hopf-Lax formula is evaluated. The Execution cost is thus 0(iVlniV). Combining the costs of 
Preprocessing and Execution, we obtain as announced the complexity (D(NlnN + Nlnn(M)). 

Remark 2.5 (Memory requirements). The memory requirements of numerical methods for 
the eikonal equation, such as the AGSI and the FM-LBR, are dominated by (I) storing the 
riemannian metric M, sampled on the discrete domain Q*, and the discrete solution d, and (II) 
storing the graph structure underlying the numerical scheme. Point (I) requires two tables of 
Nd(d + l)/2 and N reals, typically represented in 64bit floating point format, independently of 
the method. 

Point (II) is avoided for AGSI when this method is executed on a mesh with a trivial periodic 
structure, which is the case in our experiments. For the FM-LBR, Point (II) is divided into the 
storage of the direct stencils [V*(x)) X £Q n , and of the reversed stencils (V*(y)) y< zn n . Storing a 
directed graph with N vertices and N' edges requires a table of N' vector.^ and another of N 
integers: the indices of the start of the consecutive stencils in the previous table. For the direct 
or the reversed stencils of the FM-LBR, N' = 6N in two dimensions, and N' = 14N in three 
dimensions, see Propositions 1.9 and \l.l(\ We represented integers in 3 2bit format, and vector 



components in 8bit format, since these are small integers by construction. Summing up, we find 
that the memory requirements of the FM-LBR are approximately twice those of the AGSI (for 
which we do not count any mesh structure). 



The edge xx' may also be represented by the index of the endpoint x' , instead of the associated vector. 
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2.1 Convergence analysis 



We establish in this section the convergence result announced in Point (iii) of Theorem 2.2 The 
FM-LBR is applied for each n > 2 to a fixed Riemannian metric A4 G C°(Q, St), and produces 



a sequence of maps d n : Q n — > R + satisfying (10). 



Lemma 2.6. • (Limited extension) There exists a constant Vq = Vo(Ai) such that each 
vertex v of each mesh T{z), z£fi, satisfies \\v\\ < Vq. 

• (Consistency) There exists a constant ro = ro(.M) > such that the following holds. For 
each z G ^ there exists a basis (u\, ■ ■ ■ , u^) of 7L d such that u\, ■ ■ ■ , Ud, —U\, • • • , — Ud are 
vertices of the mesh T{z + w), for all w G M d such that \\w\\ < ro. 



Proof. We prove this lemma for the mesh constructions given in Propositions |1.9| and 1.10 The 
"block diagonal" case immediately follows by considering each block separately. 



The elements of a M-reduced basis of TL have a norm bounded by k(M), see Proposition 
1.6 Hence any vertex of a M-reduced mesh T built as described in Proposition |1.9| satisfies 



| < 2k{M) (resp. Proposition 1.10 and < 3k(M)). This establishes Point (Consistency). 
We denote by ro = ro (Ai ) the largest positive constant such that for all z G f2 and all w G TR d 

\w\<r => d x (M(z),M(z + w)) <ln(l + k(M)- 2 )/ 4. 



If z G f2 and \\w\\ < ro, then it follows from Corollary 1 1 . 8| that the vertices of T(z + w), which is 
Ai(z + w)-reduced mesh, contain as a subset any A1(z)-reduced basis (u\, • • • , ua) of 7L d . This 
establishes Point (Limited extension), and concludes the proof. □ 

Our next lemma shows that the discrete maps d n : Vt n —¥ 1R+ obey a Lipschitz regularity 
property, if n is sufficiently large, and thus take finite values, 

Lemma 2.7. For any n > 0, any z G £l n and any vertex v ofT{z), one has 

d n (z) - d n (z + v/n) < \\v\\ M (z)/ n - ( 34 ) 
For any n > uq, any z G Q n , and any v G { — 1, 0, l} d one has 

\d n (z) -d n (z + v/n)\ < A /n, (35) 
where n := d 2 V d /r , A := d 2 V d M , and M := max{^/\\M{z)\\; z EO,}. 



Proof. If z G O n \ {0}, then (34) follows from the equality d n {z) = A n (d n ,z), and from the 



definition (27) of A n (d n , •). If z = 0, then the left hand side of (34) is negative, while the right 



hand side is non-negative. This concludes the proof of the first part of this lemma. 



We next turn to the proof of (35 ), and for that purpose we c onsider a fixed n > uq. Consider 
a point z G f2 n and a basis (ui, • • • , Ud) of 7L d as in Lemma 2.6 (Consistency), and observe that 



ll^iH < Vo for all 1 < i < d. Let A be the d x d matrix which columns are u\, • • • , Ud- We denote 
by com(A) the comatrix of A, and we observe that the coefficients of this matrix are bounded 
in absolute value by Vg d_1 (indeed they are determinants of (d — 1) x (d — 1) sub-matrices of 
A, the norm of which columns is bounded by Vq). Since | det^4| = 1, the absolute value of the 
coefficients of A~ 1 = com(A) T / detA is also bounded by Vq~ 1 . 

Let (ai, • • • , ad) G 7L d be such that v = a±ui + • • • + ayui, in other words (a±, • • • , ay) = 
A~ 1 v. We denote s := |ai| + • • • + \ad\, and we observe that s < cPVq 1 , since the absolute 
value of the coefficients of v is bounded by 1. There exists vq,v\,--- ,v s G 7L d , such that 
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vq = 0, v s = v and Vi+% — V{ G {«!,-•• , Wd, — ui,--- , — lid} for all 1 < i < s. Furthermore 
\\vi\\ < sVo < c^ 2 ^' hence ui, • • • , u d , — «i, • • • , — u d are vertices of Tiz + Vi/n) according to 
Lemma |2.6| (Consistency). It thus follows from (34) that 



n(d n {z + Vi/n) - d n (z + v i+1 /n)) < \\v i+1 - Vi\\ M ( z+Vi ) < \\v i+1 - Vi\\y/\\M(z + Vi)\\ < M Q V , 

which implies that d n (z) — d n (z + v/n) < sMqVo/ti < cPMqVq /n = A$/n. Exchanging the roles 
of z and z + v/n we obtain likewise d n {z-\-v/n) — d n {z) < Aq/u, which concludes the proof. □ 

Let ip : lR d — > f2 be the canonical surjection. We denote by d per the distance on defined 
for all z, z' G by 

d P er(^, z) := min{||Z - Z'\\ Z £ f^ l {z), Z' G ^(z 1 )} 
Corollary 2.8. For each n > no and for all z GO we define 

D n (x) := min d n (z) + A d per (z, x), (36) 



where the constants no, Ao, are defined in Lemma 2.1 

The map D n : (SI, d per ) — > JR+ is Ao-Lipschitz, and d n (x) = D n (x) for all x G tt n . 

Proof. The map D n is Ao-Lipschitz since it is defined as the minimum of the family of Ao- 
Lipschitz functions x \-t d n (z) + Ao d pcr (z, x), indexed by x G Q n . 

Let x G O n , let v = (fi, • • • ,v d ) G 7L d and let v mSiX := max{|ui|, • • • , \v d \}. Applying v max 
times (35) we obtain 

| d n (x) - d n (x + v/n)\ < Ao^max/"- < A ||f/n||. 

It follows that d n (x) < d n (z) + Aod per (x,z) for all x,z G O n . This immediately implies as 
announced that d n (x) = D n (x) for all x G O n , which concludes the proof. □ 

Before stating the main result of this section, we recall the definition of the viscosity solution 
of an eikonal equation |llj . 

Definition 2.9. The viscosity solution of the eikonal equation Q), is the unique continuous 
function D : O — > M + such that for any ip G C 1 (0, M) and any z G O \ {0} the following holds: 

• IfD—ip has a unique global maximum at z, then ||Vy(2)||^vt( z )-i < 1. 

• IfD—{p has a unique global minimum at z, then || V</?(^)||x(z)- 1 — 1- 

In alternative definitions of the notion of viscosity solution, the assumption that "D —ip has 
a unique global maximum at z" is often replaced with "D —ip attains a local maximum at z" 
(resp. minimum). These two definitions are equivalent, since one may subtract (resp. add) to 
if a suitable smooth function ip G C 1 (0,]R + ), large far from z, and with a parabolic behavior 
close to z: ip(z + h) ~ A||/i|| 2 for h sufficiently small. 



We show in the following that the functions D n , introduced in Lemma 2.8 converge uniformly 
towards the viscosity solution of ^ as n — > oo. This immediately implies the convergence of 
the discrete maps d n produced by our modified algorithm, as stated in Point (iii) of Theorem 



2.2 The proof is similar in essence to the proof provided in [5], yet with a number of minor 



modifications due to our specific context. 
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Since (f2,d per ) is a compact metric space, since D n (0) = and since D n is Ao-Lipschitz for 
all n > uq, the theorem of Arzela-Ascoli implies that the sequence D n is pre-compact. In the 
rest of the proof we consider an arbitrary subsequence (D ff ( n ))„>o which converges uniformly 
on £1 towards a limit D : — > 1R + . Our objective is to establish that D is a viscosity solution 
of the eikonal equation ([2]). Once this point is established, we conclude as announced that D n 
converges uniformly toward the viscosity solution of Q, using the uniqueness of the viscosity 
solution and the pre-compactness of the sequence D n . 

We consider for each n > 1 a point z n G Q n , the non-zero vertices v™, • • • , v 1 ) G TL d of a 
common simplex T G T(z n ), and some coefficients a", • • • , G H + , a" + • • • + = 1, which 
will all be specified later. We consider A n > 0, and v n G M, d , \\v n \\ = 1, such that 

X n v n = « + • • • + a n d v n d . 



Note that X n < Vq. On the other hand using the acuteness property (13) and denoting M :- 
M.{z n ) we obtain 

Anl|M|| > ||A^n|||f 

l<i<n l<*<i<^ 

min \\Vi\\ 2 M 

Ki<n 




hence A n > {cLk{M.)) 2. Denoting 

D n := d n (z n ) - a? d n {z n + v?/n) 

Ki<d 



it follows from the definition (27) of the Hopf-Lax update that 

D n < X n \\v n \\M{z n )/ n (37) 

but also that, given n > and z n E we may choose v • • • , v d and a™ , • • • , a d m sucn wa Y 
that the above inequality is an equality. 

Consider an arbitrary map ip G C 1 (17,11) and a point z G fi \ {0}. Denoting L := V</?(,z) T 
we have for any x G f2 and any /i G M d the Taylor development 

\tp(x + h) - (p(x) - Lh\ < [ \\V<p{x + th)- V(p(z)\\\\h\\dt<u(d p0r (x,z) + \\h\\)\\h\\, (38) 

Jo 

where uj denotes the modulus of continuity of the continuous function Vy> : f2 — > TR d . We denote 

S n := (p(z n ) - afyC^n + ^ n / re ) 5 

l<i<d 

and we observe that, using ( |3~8~| ) and denoting r n := d per (^, z n ) + Vb/n, 

|5„ + \ n Lv n /n\ < Y a i Wi z n + v?/n) - <p(z n ) - Lv?/n\ < ui(r n )V /n. (39) 

KKn 
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We now explicit the choice of z n , (vf) and (af ), and we use the above inequalities to conclude 
the proof. Let us assume that D — ip has a strict global maximum (resp. minimum) at z. For 
each ii > we denote by z n G £l n a maximizer (resp. minimizer) of D — d n on Q n . It follows 
from the uniform convergence T) a r n \ — >• D, as n — > oo, that z a r n \ — > z. 

We consider a fixed v 6 M d such that ||v|| = 1, and we choose a", • • • , and vf, ■ ■ ■ , as 
above in such way that v n = v for all n > 0. (resp. We choose a", • • • , and t>™, • • • , Vj in 
such way that ( |37[ ) is an equality for all n > 0, and we denote by v an arbitrary cluster value of 
the sequence (u CT ( n ))n>2-) 

Using successively (39), the definition of z n , and (37) we obtain 

-\ n Lv n /n - uj(r n )V /n < S n < D n < X n \\v n \\M(z n )/ n 

(resp. -X n Lv n /n + u(r n )V /n > S n > D n = K\\vn\\M(z n )/n)- I* follows that 

-Lv n < \\v n \\ M ( Zn ) + u(r n )V /\ n 

(resp. -Lv n > ||u n ||.M(z„) ~ u(r n )Vo/\ n ). Since r CT ( n) as n -)■ oo (resp. r CT ( n ) -)■ and u is 
a cluster value of the sequence (v a ( n )) n> o), and since A n is bounded below independently of n, 
we obtain 

-Lv < \\v\\ M{z) 

(resp. — Lv > \\v \\m{z))- Observing that 

l|Vy(z)|U(*)-i =sup| ^~^ ; u G ]R rf , ||v|| = l| 

we obtain as announced that ||V(^(z)||_yv(( 2 ) < 1 (resp. ||V(/?(2;)||_A4( 2 ) > 1). It follows that D is 
the viscosity solution of Q, which concludes the proof. 



3 Numerical experiments 

This section is devoted to the numerical comparison of three solvers of the eikonal equation: 
two popular methods (AGSI, FM-8) which enjoy a reputation of simplicity and efficiency in 
applications, and the proposed algorithm (FM-LBR). The Adaptive Gauss Seidel Iteration^] 
(AGSI) [5] produces numerical approximations which are guaranteed to converge towards the 
solution of the continuous anisotropic eikonal equation as one refines the computation gricQ for 
an arbitrary continuous Riemannian metric A4. Fast Marching using the 8 point stencil (FM-8, 
stencil illustrated on Figure [TJ center left) does not offer this convergence guarantee, but has 
a quasi-linear complexity 0(N In N), in contrast to the super-linear complexity 0{fi{Ai)N 1+ d^ 
of the AGSI. Fast Marching using Lattice Basis Reduction^] (FM-LBR) aims to offer the best of 
both worlds: a convergence guarantee, and fast computation time^J 

2 As suggested in [5], the stopping criterion tolerance for the iterations of the AGSI is set to 10~ 8 . 

3 The grid is triangulated with a trivial periodic mesh, for the AGSI and the MAOUM 

4 



We use for each discrete point z the A4(«)-reduced neighborhood described by Proposition 



1.9 



in 2-d, and 

Proposition 1.10 in 3-d, except if the matrix A4(z) is detected to be exactly diagonal. In that case we use the 
standard 4 vertices neighborhood in 2-d (resp. 6 vertices in 3-d), which is a M (z)-reduced mesh, see Figure [l](left 
and center right). This modification has little impact on accuracy or CPU time, but avoids to pointlessly break 
the symmetry of the numerical scheme. A C++ source code, provided as an ancillary file to the Arxiv version of 
this paper, allows to reproduce the above experiments. 

5 Note that memory requirements are doubled for the FM-LBR in comparison with the AGSI and FM-8, see 



Remark 2.5 
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First test First test, rotated by x/6 Second test Third test 




-0.4 -0.2 0.0 0.2 0.4 -0.4 -0.2 0.0 0.2 0.4 -0.4 -0.2 0.0 0.2 0.4 -0.4 -0.2 0.0 0.2 0.4 



Figure 6: Level lines of the solutions of the two dimensional test cases. 



We also implemented the Monotone Acceptance Ordered Upwind Method (MAOUM) [2], a 
Dijkstra inspired method using static stencils, like the FM-LBR. The difference between these 
two methods is that the stencil V n (z) used by the MAOUM at a point grid point z £ Q n 
is isotropic, only depends on the anisotropy ratio k(A4(z)), and its boundary has cardinality 
#(f2 n n dV n {z)) = 0(k(M(z))); in contrast the stencil of the FM-LBR is anisotropic, aligned 
with the ellipse defined by M(z), and of cardinality 0(1). The stencils of the MAOUM were 
precomputed and stored in a look-up table, resulting in a complexity 0(n(M)N In N) for this 
algorithm. 

We consider four test cases. The first two involve two dimensional Riemannian metrics A4 
of limited anisotropy, k(M) ^ 5, and were proposed in [231 [19] (including the grid sizes). All 
four methods (FM-LBR, FM-8, AGSI, MAOUM) are viable, yet the proposed method shows 
a reduced CPU usage and a competitive accuracy in comparison with its alternatives. The 
third test case, introduced in [3] and related to image segmentation problems, involves a highly 
anisotropic two dimensional metric M, k(M) ~ 100. It was observed in [3] that neither the 
AGSI or the FM-8 are viable in this context, due to complexity and accuracy issues respectively; 
and the MAOUM is no better. We show that FM-LBR solves this dilemma. The fourth test is 
a generalization of the third one in three dimensions. 

Remark 3.1. Contrary to the assumptions of the convergence analysis, in but consistently 
with the envisioned applications, our test cases do not involve periodic boundary conditions. As 
a result, in exceptional circumstances, a few grid points usually in the domain's corners may not 
be reached by the FM-LBR, because they are not connected to the origin in the underlying graph 



(Lemma 2.1 does not apply). In the following experiments this only happened when the first test 
case was rotated by an angle 9 £ [0.56,0.61] radians. The numerical values at the four corners 
of the grid, equal to +oo, were then rejected when computing the L 1 and L°° errors. 

The last two test cases involve discontinuous Riemannian metrics, which happens in many 
applications (e.g. at the junction between materials of different index in geometrical optics), but 
also contradicts the assumptions of our convergence analysis. 

The first benchmark, introduced in |24| , consists in computing the Riemannian distance from 
the origin (0,0,0) on the surface parametrized by / : (x,y) — > (x, y, (3/4) sin(37rx) sin(37ry)). 
Denoting by F := (d x f, d y f) the differential of /, the Riemannian metric is given by M. = 
F^F. The coordinates (x,y) are restricted to the square [— 0.5,0. 5] 2 , and the anisotropy ratio 
is k(A4) ~ 5.1 [23]. A reference solution is obtained on a 4000 x 4000 grid using the AGSI, and 
extended by bilinear interpolation. The FM-LBR, FM-8, AGSI and MAOUM are applied on 
a 292 x 292 grid, and the L°° and averaged L 1 errors estimated with respect to the reference 
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Figure 7: Tables of CPU time in seconds, L°° error an averaged L 1 error (left). Accuracy, in 
the first test rotated by an angle 9 G [0, 7r/4] (this interval is enough, since the dependence in 
9 is 7r/2-periodic and even). In average over theta, CPU times are 0.21s, 0.20s, 1.37s, 1.31s, 
L°° errors 5.16, 7.64, 6.86, 8.57 and averaged L 1 errors 1.34, 2.58, 1.95, 2.40 for the FM-LBR, 
FM-8, AGSI and MAOUM respectively. All errors are multiplied by 100, for better readability. 

solution. As shown on Figure [7J FM-LBR is the fastest in terms of CPU tim^J but is less 
accurate than the AGSI or the FM-8. 

Rotating the this test case by the angle 9 = vr/6, and conducting the same experiment, shows 
a different story: the numerical error increases strongly for the AGSI and the FM-8, while the 
FM-LBR, unaffected, is now the most accurate method, see Figure [7| Unsurprisingly the CPU 
time is almost unaffected for the one pass solvers (FM-LBR, FM-8, MAOUM), while it slightly 
increases for the AGSI. The FM-LBR cuts L°° and L 1 numerical errors by 40% in comparison 
with the AGSI and the MAOUM, and CPU time by 85%, while the FM-8 produces even larger 
numerical errors. 

Figure [7] shows that the FM-LBR offers the best accuracy for more grid orientations 9 than 
its alternatives. The maximal error and averaged error with respect to 9 are also in favor of the 
FM-LBR. The heuristic explanation of these empirical observations is the following: for 9 = 0, 
this test case is dominated by (close to) axis-aligned anisotropy. The fixed stencils of the AGSI 
and the FM-8 seem to benefit from this configuration; the FM-8 also works well for 9 = 7r/4, be- 
cause its stencil includes the four diagonals. In contrast, the large stencils used by the MAOUM 
provide a good angular resolution in all directions, hence the method is mostly unaffected by 
rotations of the discretization grid. The FM-LBR does not suffer either from off-axis anisotropy, 
thanks to its adaptive stencils which are aligned with the anisotropy directions. In the special 
case of constant Riemannian metrics, random grid orientations are actually more favorable in 
average to the FM-LBR than axis aligned ones, see |12j . 

The second benchmark, discussed in [24|, I19j. is inspired by seismic imaging. The Rieman- 
nian metric M is defined as follows: at each point z = (x,y), the symmetric matrix A4(z) has 

6 All timings obtained on a 2.4Ghz Core 2 Duo, using a single core. Timings of the FM-LBR include the stencil 
construction, which typically accounts for 25%. 
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Figure 8: Reference solution for the third test case (left). The Riemannian metric M is 
anisotropic only on a thin band along a spiraling curve, wide of a few grid points. Detail at 
resolutions n x n, where n equals 200 (center left), 500 (center right) and 1000 (right). 



the two eigenvalues 0.8, 0.2, the latter associated to the eigenvector (1, (n/2) cos(4-7rx)). There 
is little or no bias here towards axis-aligned anisotropy. The Riemannian distance from the 
origin (0,0) is computed on the square domain [0.5, 0.5] 2 . A reference solution is obtained on a 
4000 x 4000 grid using the AGSI, and compared with the results of the four algorithms, obtained 
on a 193 x 193 grid. As shown on the table Figure [7[ FM-LBR takes a smaller CPU time and 
offers a better accuracy than its alternatives. 



The third [I] and fourth test cases are relevant benchmarks if one's objective is to use fast 
marching methods for the segmentation of tubular structures, in medical image or volume data 
respectively. The Riemannian metric M is euclidean (equal to Id) on a box domain Q, except 
on the neighborhood of a curve T, see Figures [8] and 12 (right). The matrix A4(z), when z € O, 



is close to r, has two eigenspaces: one of dimension 1, directed "tangentially" to the curve T, 
and associated to a small eigenvalue (1/100 2 and 1/50 2 respectively), and one of co-dimension 
1, associated to the eigenvalue 1. The (approximated) Riemannian distance D to to the origin 
G is computed, and a path of minimal minimal length 7 joining a given point P 6 is 
extracted by "gradient descent on the Riemannian manifold (£l,A4)": 

7 , (t) = -A4( 7 (t))- 1 VD( 7 (t)). (40) 

By construction of the Riemannian metric traveling close and tangentially to the curve V is 
cheap. This is reflected by the level lines of D, and by the allure of the minimal path, see Figures 
[6j|8|and|12| Heuristically, this path joins the curve T in straight line, almost orthogonally, and 
then follows it. The alignment of the minimal path with the direction of anisotropy, observed in 
this test case, is not an uncommon phenomenon. The FM-LBR presumably benefits a lot from 
this behavior in terms of accuracy, since its stencils typically provide a good "angular resolution" 
in the direction of anisotropy, see Figures [2j [3j |5[ 

The performance of the FM-LBR, FM-8 and AGSI is illustrated on Figure [9] for the third 
test case. The MAOUM showed a poor accuracy in this test, presumably due to the huge sten- 
cils it generated, hence its results are not shown. The CPU time/resolution curve of the AGSI 
shows a stronger slope than the one pass solvers (FM-LBR, FM-8), which reflects its intrinsically 
larger complexity, namely 0(n(Ai)N 1+ ~3^ instead of 0(N In N). The L°° and L 1 error curves 
suggest that the FM-8 is not consistent in this test case, contrary to the FM-LBR and AGSI. 
The reference solution was obtained on a 4000 x 4000 grid, using the variant of the Fast March- 
ing algorithm described in [13]. The AGSI could not be used due to prohibitive computation 
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Figure 9: CPU Time (left, in seconds), L°° error (center), and averaged L 1 error (right) of the 
FM-LBR, FM-8 and AGSI, at several resolutions ranging from 120 to 1200 (log-log scale). 

times and unimpressive accuracy in this test, and the FM-LBR was rejected because using it 
for reference might have induced a bias in our favor. At the resolution 1000 x 1000, typical in 
image analysis, the FM-LBR cuts the L°° error by 80% and the L 1 error by 75% with respect 
to the AGSI, while reducing the CPU time from 11 minutes to 2.5 seconds. 

Obtaining the shortest path joining two given points is essential in motion planning control 
problems [2], as well as in the envisioned application [4J. This involves solving the Ordinary 



Differential Equation (ODE) (40), a task less trivial than it seems. The author is conscious that 
a benchmark based on minimal paths may reflect the properties of the ODE solver (and the time 
spent adjusting its sometimes numerous parameters), as much as those of the eikonal solver, but 
does so nevertheless due to the importance of minimal paths in applications. Eikonal solvers 



based on the discrete fixed point problem (10), such as the FM-LBR, FM-8 and AGSI, provide 
at each grid point z G £l n an estimate d(z) of the distance D(z), and in addition an estimate 
v(z) of the direction and orientation of the distorted gradient —M.(z)~ 1 VD(z). This estimate 
has the form ^ 

v ( z ) '■= - ^2 aiVi = a ^ Zi ~ ( 41 ) 

l<i<k l<i<k 

where the integer 1 < k < d, the positive barycentric coefficients (Qj)^ =1 and the vertices (v i)f =1 
of the mesh T(z) (resp. vertices {zi)\ =l of V n (z)) achieve the minimum in the Hopf-Lax update 



operator (27). 



From this point, a typical approach to solve (40) is to extend the values of d or v to the 
continuous domain Q via an interpolation procedure, and then to use a black box ODE solver or 
a Runge Kutta method. Note that the accuracy usually expected from these high order methods 



is mostly doomed, since the discretization ( 10 ) of the eikonal equation is only first order, and 



since the vector field M _1 VD is discontinuous both at "caustics" and discontinuities of M. A 
more significant issue is that computations frequently get stuck, despite the use of state of the 
art and/or commercial interpolation methods and ODE solvers, see e.g. [2] Figure 5.10. 

We propose a method for the computation of minimal paths, which trades high order accuracy 
for robustness, and never gets stuck if the eikonal solver is Dijkstra inspired. It takes advantage 



of the specific form ( 10 ) of the discretization of the eikonal equation, and does not rely on black 



box routines. It is parameter free: there is not interpolation order or gradient step to adjust. 



As illustrated on Figure 1 1 , the FM-LBR recovers a qualitative minimal path in the third 
test case for grid resolutions much coarser than the AGSI. See Figure 8 for an illustration of 
the impact of coarse resolutions on the problem discretization. Hence the better accuracy of 
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Figure 10: Notations for the minimal path computation (left); the contour of the stencil V n (xi) 
is shown dotted. Grid points xq, ■ ■ ■ , x r G £l n , corrections uq, ■ ■ ■ ,u r G H d shown as arrows, and 
piecewise linear path 7 (center). Grid Points (xj)[ =1 and piecewise linear path 7 in the second 
test case at resolution n x n, n = 100, using the FM-LBR (center right, detail). 

the FM-LBR in this test case, see Figure [9j does effectively translate into a better extraction of 
minimal paths. 

Our method produces a piecewise linear path 7, joining the points xq + uq, • • • ,x r + u r of 
fi, where Xq, ••• ,x r G £l n are grid points, and Uq, ■ ■ ■ ,u r G lR d are small correcting offsets, 



see Figure 10. The first grid point xq is given by the user (and should satisfy d(xo) < 00, see 
Remark 3.1 ), and the first offset u$ is zero. All produced grid points Xi+i, < i < r, are vertices 
of the stencil V n {xi) of the previous point. In the case of a Fast Marching method, the causality 
property implies that the values (d(xj))[ =0 are strictly decreasing, down to zero: d(x r ) = (the 
last point belongs to the source). 

Minimal path computation, starting from a given grid point P G Q n . 
Initialisation: xq := P, uq := 0, i ::= 0. 
While d(zi) > do 



Denote by z\, ■ ■ ■ ,Zk the grid points appearing in the expression (41 ) of v(xj). 
Find A G H+, 1 < j < k, which minimize \\xi + u% + Xv(xi) — Zj\\ 
Set Xi + i := Zj and Ui + \ := xi + Ui + Xv{xi) — Zj. 
i ::=i + 1. 



The piecewise linear path 7 : [0, r] 
differential equation 



fi, parametrized so that 7(2) 
j'(t) = X [tl v(x [tl ), 



+ Ui, satisfies the 



(42) 



for non-integer t G [0,r], where the constants (Aj)|" =1 are the minimizers in the second step of 
the while loop, and the vector v(z), z G Q n , is defined in (41). The particularity of our path 
extraction method is that the direction field v is not evaluated on the curve 7, but at the nearby 
points (xj)£_0. The next proposition shows that these points remain at distance 0(k(A4) 3 /n) 
from the curve in dimension 2 (the exponent 3 seems to be either over-estimated or a rare 
worst case scenario, in view of the experiments, see Figure 10), hence the accuracy of our path 
extraction method should be no worse than an explicit Euler integration of the ODE (40), with 
step n(M.) 3 /n. We only use the fact that the meshes T(z), zfll, have diameter 0(k(A4)), and 
are built of triangles of area 1/2. 

Proposition 3.2. Consider a piecewise linear path 7 G C°([0,r],Sl) extracted with the above 
algorithm and parametrized as in (4%), in dimension 2, after the Execution of the FM-LBR on 
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Q n for a Riemannian metric A4. Then for all t € [0, r\, one has for some absolute constant C 

h(t)-x [t] \\<CK(Mf/n. (43) 

Proof. We consider a fixed integer < i < r, and observe that for all t E [i, i + 1[ one has, since 
7 is linear on this interval and 7(2) = Xj + Ui 

||7(t) - x\t] II < max{||7(i + 1) - Xj||, ||7(i) - Xj||} < max{||u i+ i|| + \\x i+ i - Xi\\, \\ui\\}. 



By construction, there exists a non-zero vertex Vj of the mesh T(xi), appearing in (41), such 



that Xi+i — Xi = Vj/n. Combining Propositions 1.6 and 1.9 we obtain that ||u,-|| < 2k(A4), hence 
\\xi+i — Xi\\ < 2k(M)/u. 

Our next objective is to bound uniformly. In order to avoid notational clutter, we 
rescale the variables of interest by a factor n: we set u := niti, u' := m/j+i, and v := nv(xi). By 
construction (second step of the while loop in the path computation), we have for any A G H+ 
and any 1 < j < k 

\\u'\\ < \\u + Xv - Vj\\, (44) 
where the integer 1 < k < 2, and the vectors (1^)^=1 are those appearing in the expression 



(41) of v. If k = 1, then choosing A = 1, j = 1, and observing that v = v±, we obtain 
||w'|| < 1 1 1* + 1 X V — Vl\\ = \\Ui\\. 

We next assume that k = 2, and observe that v = a\V\ + 0^2 for some «i,«2 S 1R+ such 



that «i + 02 = 1, see (41). We consider an arbitrary but fixed [i G]l,oo[ (e.g. fi := \/2), and 
define w\ := ui// — ^2//", u>2 := — Note that \\wj\\ < (/x + /x -1 ) max{||ui||, ||f 2 1| } — 

Co«(-M) with Co = 2(/i + /x _1 ). The matrix A of lines uii and W2 satisfies |detA| = (/j, 2 — 
H~ 2 )\ det(ui, ^2)! = fJp — fJ,~ 2 , since | det(ui, V2)\ = 1, see Point (b) of Definition |l.l[ On the 
other hand \\A\\ < Vlkill 2 + IMP < \/2C k(M), thus [|^4_— 1 [|— 1 = | det A|/||A|| > C x /k(M) 
where C\ := (/i — /.i _1 )/(v / 2Co). Therefore there exists 1 < j < 2 and e E {—1, 1} such that 



T 

eu - 



: wjV2 > ^(u^wj) 2 + (u T Wj ) 2 = \\Au\\ > H^ir^MI > Cx\\u\\/k(M). 
Assume, without loss of generality, that j = 1 and e = 1 satisfy the above expression. Then 



choosing A = l/(«i + \x a.2) and j = 1 in (44) yields, with v := 02 A 



" |2 < ||u - v Wl \\ 2 = llti.il 2 - 2^ T wi + 4i/ 2 ||u;i|| 2 < ||ti|| 2 - 2vCi\\u\\/k(M) + v 2 C$k(M) 2 



If ||it|| > C2k(M) 3 , where Ci ■= Cq/(2Ci), then the above inequality shows that < 



since v < 1. If is below this bound, then choosing A = in (44) yields \v!\ < \\u\\ + < 



C2k(A4) + 2k(A4). By an immediate induction argument we thus obtain ||«j|| < (C2k(M) + 
2k(M))/ti, for all < i < r, which concludes the proof. □ 

The fourth test case is a generalization of the third one to three dimensions. It also involves 
a spiraling curve T, surrounded by a thin tubular neighborhood where traveling tangentially 
to r is 50 times cheaper than in orthogonal directions, or anywhere else in the domain, see 



Remark 3.3 for details. CPU time was 105s for the FM-LBR, while the AGSI took 480s and 



failed to recover the minimal path presented on Figure 12 (center) (a straight line joining the 
two endpoints was obtained instead). We do not perform a detailed benchmark in this case, 
which is similar in essence to the third one. Let us simply point out that FM-LBR addresses 
here a large scale (more than 10 millions grid points), strongly anisotropic (k(A4) = 50) three 
dimensional shortest path problem, with a good accuracy and within reasonable CPU time on 
standard laptop computer. 
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Figure 11: Visual comparison of the accuracy of three algorithms, at three resolutions, in the 
2-d test case. Qualitatively, the approximate geodesic has the right behavior for a resolution as 
low as 170 x 170 with the FM-LBR, and 1000 x 1000 with the AGSI. This is presumably never 
the case for the FM-8, which is not consistent here. 




Figure 12: Results of the FM-LBR in the fourth, 3-d, test case. Iso-surface {d(z) = 2} (left), 
and shortest path joining the points (0,0,0) and (3,0,0) (center). Detail of the discrete points 
(represented by small cubes), in the neighborhood of the curve T(t) = (coswoi, sinwoi, t), for 
which the Riemannian metric is not euclidean (right). 
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Remark 3.3. Exact description of third and fourth test cases, given for reference. 

• The third, 2d test case is posed on the domain Q = [—1.1, 1.1] 2 , involves three parameters 
ujq = 6ir, ro = <5o = 1/100 (the parameter 5o was it seems slightly larger in the original 
experiment ^j), and the curve F : [0, 1] — > O defined by T(t) = t(coscoot, sinwoi). For all 
z G £1 we have Ai(z) = Id, except if there exists < t < 1 and < r < r^ such that 
z = T(t) + r(cosujot, sinwo^)- In that case Ai(z) has the eigenvalues 5q and 1, the former 
being associated to the eigenvector T'(t). Geodesic starting point: P = (1, —1). 

• The fourth, 3d test case is posed on the domain Q = [— 1.1,1. 1] 2 x [0,3], discretized on 
a 200 x 200 x 272 grid, involves the three parameters ujq = (5/2)ir, ro = 5$ = 1/50, and 
the curve defined by T(t) = (cos uot, sin uot, t). For all z £ Q we have M.(z) = Id, except 
if there exists t,\,fi £ M such that z = F(t) + (A cos(ujot), X sinwot, //) and A 2 + /x 2 < 
(r /2) 2 . In that case Ai(z) has the eigenvalues 8q and 1, the former being associated 
to the eigenvector T'(t) and the latter to the orthogonal space. Geodesic starting point: 
P = (0,0,3). 

Conclusion 

In this paper, we introduced a new discretization scheme for anisotropic eikonal equations: Fast 
Marching using Lattice Basis Reduction (FM-LBR) . It is a variant of the classical Fast Marching 
algorithm, based on the algebraic concept of Lattice Basis Reduction, which strongpoints are 
the following. (I, Convergence) This algorithm is consistent for the anisotropic eikonal equation 
associated to any continuous Riemannian metric, of arbitrary anisotropy. (II, Complexity) It has 
a numerical cost comparable to classical isotropic Fast Marching, independently of the problem 
anisotropy. (Ill, Accuracy) The accuracy of the FM-LBR is competitive in general, and striking 
in test cases, related to tubular segmentation in medical images, where the Riemannian metric 
has a pronounced anisotropy close to and tangentially to a curve. 

These strongpoints come at the price of the specialization of the FM-LBR: (i) the Riemannian 
metric may not be replaced with a more general Finsler metric, (ii) the domain needs to be 
discretized on a cartesian grid, and (iii) of dimension 2, 3, or in higher dimension the underlying 
Riemannian metric needs to have a block diagonal structure. Hopefully these requirements are 
met in many applications, and future work will be devoted to the application of the proposed 
algorithm in the context of medical image processing. 
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